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THE  EFFECT  OF  ATMOSPHERIC  TURBULENCE  ON  AIR  TO  GROUND  PHOTOGRAPHY 

by 
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SUMMARY 

The  current  state  of  the  theory  of  optical  effects  of  atmospheric  turbu¬ 
lence  is  summarized,  and  the  derivation  of  an  <  tuation  relating  the  modulation 
transfer  function  of  an  optical  path  through  the  atmosphere  to  a  weighted  inte¬ 
gral  of  the  structure  constant  of  refractive  index  along  that  path  is  outlined. 

Methods  of  estimating  the  structure  constant  from  more  readily  available 
meteorological  parameters  are  discussed,  and  FORTRAN  subroutines  for  evaluating 
the  integral  and  estimating  the  modulation  transfer  function  of  the  path  under 
any  required  conditions  are  given. 
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1  INTRODUCTION 

When  the  optical  path  of  an  image  forming  system  passes  through  a  fluid  medium*  the 
random  variations  in  its  refractive  index  arising  from  turbulence  modify  the  image  formed 
in  three  main  ways. 

(1)  The  total  light  flux  in  the  image  of  a  point  source  varies  with  time. 

(2)  The  image  of  the  point  is  spread  out  into  a  patch  of  light. 

(3)  This  patch  moves  about  randomly  in  the  image  plane. 

The  first  of  these  effects  produces  the  scintillation  of  the  stars*  but  is  not  usually  of 
importance  when  imaging  extended  sources.  The  second  is  the  most  important  when  photo¬ 
graphing  terrestrial  objects  with  short  exposure  times*  but  with  longer  exposures*  the 
motion  of  the  aerial  image  also  contributes  to  the  spread  of  the  recorded  image. 

In  aerial  photography  two  types  of  turbulence  are  of  importance*  the  local  effects 
produced  by  the  aircraft  itself  as  it  passes  through  the  air*  and  those  occurring  through¬ 
out  the  whole  of  the  optical  path  from  meteorological  causes.  Either  of  these  can  limit 
the  performance  of  an  imaging  system,  but  the  former  is  highly  dependent  on  the  precise 
location  of  the  camera  in  the  aircraft*  and  on  its  local  and  overall  aerodynamic  design. 
Every  case  must  be  considered  individually  and  experimental  measurement  on  an  actual 
installation  is  probably  the  only  satisfactory  way  of  determining  the  magnitude  of  the 
effect.  In  this  Report  we  will  be  concerned  solely  with  the  meteorological  effects. 

Two  distinct  types  of  problem  are  involved,  estimating  the  variations  of  refractive 
index  along  the  optical  path*  and  computing  the  effect  of  these  on  the  quality  of  the 
image  formed  by  an  optical  system  at  the  end  of  that  path.  Both  of  these  present  formid¬ 
able  difficulties. 

Experimental  measurements  of  the  turbulence  in  the  free  atmosphere  have  been 
♦  1-9 

published  by  a  numer  of  workers  #  These  show  that  the  magnitude  of  the  effect  can  vary 
very  rapidly  in  both  space  and  time*  and  it  seems  unlikely  that  it  will  ever  be  possible 
to  make  accurate  predictions  for  a  given  place  and  instant.  For  the  purposes  of  system 
performance  estimation  and  optimization*  however*  all  that  are  required  are  estimated 
rms  values  for  a  given  general  area  and  season.  It  may  ultimately  be  possible  to  base 
such  estimates  on  statistical  data  on  actual  turbulence  measurements*  but  until  a 
sufficiently  large  data  bank  has  been  accumulated  for  the  purpose*  the  best  that  can  be 
done  is  to  base  the  estimates  on  those  meteorological  parameters  for  which  systematic 
records  are  generally  available. 

The  optical  effects  of  atmospheric  turbulence  are  still  the  subject  of  intensive 
theoretical  study.  The  classical  approach  of  Tatarski^  involves  a  number  of  assumptions 
and  approximations  whose  validity  has  been  discussed*  inter  alia  by  Strohbehn'*,  but  as 
will  be  shown*  some  of  these  limitations  are  not  too  serious  in  our  particular  applica¬ 
tion,  and  any  bias  introduced  by  the  others  is  probably  small  compared  with  the  errors 
arising  from  uncertainties  in  the  refractive  index  data.  We  have  therefore  restricted 
our  treatment  to  the  classical  theory*  and  an  outline  of  the  derivation  of  an  expression 
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for  the  modulation  transfer  function  of  a  turbulent  optical  path  will  be  given  in  the 
following  section. 

2  A  SUMMARY  OF  THE  THEORY 
2. 1  Basic  concepts 

The  optical  effects  of  turbulence  arise  from  the  random  variations  in  space  and 
time  of  the  refractive  index  of  the  air  or  other  medium  through  which  the  light  passes  on 
its  way  from  the  object  to  the  sensor.  An  excellent  and  up  to  date  introduction  to  the 
subject  has  been  given  by  Hufnagel^,  and  the  treatment  given  here  will  be  based  closely 
on  this. 

The  variation  in  space  of  a  random  variable  R  may  be  described  statistically  by  a 
structure  function 


DR(r)  -  ([><£,)  -  R(E2)]2)  0) 

where  r^  and  r^  are  vectors  defining  two  points  in  space, 
and  r  »  -\  ~  -2  ' 

The  angle  brackets  denote  the  ensemble  average. 

is  related  to  the  more  familiar  parameters  oR  ,  the  variance  and  ,  the 
auto  correlation  function,  by  the  equation 

‘  2°rI  '  ”  PR(-}  f 


In  the  Kolmogerov  model  of  turbulence  in  an  isotropic  medium  the  structure  function 
takes  the  form 


dr(D 


2  * 

V3 


(2) 


.  2 

The  constant  of  proportionality  CR  is  known  as  the  structure  constant  of  R  .  It 

depends  on  the  local  properties  of  the  atmosphere,  and  may  vary  along  the  optical  path. 

Clearly  the  model  must  break  down  for  very  large  r  as,  in  a  macroscopically  homogeneous 

medium,  D  will  tend  to  a  limit  (twice  the  variance)  as  r  increases.  We  will  return 
R 

to  this  point  later. 


The  Kolmogorov  model  can  also  be  expressed  in  terms  of  the  spectral  density  function 
of  R  . 


yto 
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where  K  is  the  spatial  frequency. 


(3) 


We  will  also  make  use  of  F  the  two-dimensional  spectral  density  of  the  structure 

K 

function  of  R  which  is  related  to  it  by  the  integral  transform 
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where  p  and  k  are  position  and  frequency  vectors  in  a  plane. 

The  central  problem  is  to  relate  the  structure  function  of  the  complex  amplitude 
of  a  wave  front  at  the  aperture  of  the  detecting  system,  to  that  of  the  refractive  index 
of  the  air  or  other  medium  along  the  path  from  the  object  to  the  detector. 

The  method  used  was  developed  by  Tatarski,  who  has  given  a  very  full  treatment  of 

the  plane  wave  case*.  In  normal  terrestrial  photography,  however,  we  are  concerned  with 

spherical  wave  fronts,  but  these  can  be  treated  in  a  very  similar  way,  as  indicated, 

,  1 2 

for  example,  by  Carlson  and  Ishimaru 
2.2  The  wave  equation 

The  propagation  of  light  through  a  turbulent  medium  can  be  expressed  as  a  set  of 
scalar  wave  equations  of  the  form 

V^u  +  k^n^(r)u  «  0  (5) 

where  u  is  any  component  of  the  electromagnetic  field 
k  is  the  wave  number 

n(r)  is  the  refractive  index  at  the  point  in  space  represented  by  the  vector  r  , 

It  is  convenient  to  express  u  in  terms  of  its  log  amplitude  and  phase,  and  to 
express  these  as  the  sum  of  a  solution  to  the  unperturbed  equation  and  small  perturba¬ 
tions  to  the  solution  arising  from  the  random  departures  of  n  from  its  mean  value 
(assumed  unity) . 


u 

m 

A 

Ae 

(6) 

* 

- 

In  u  ■=  In  A  +  iS 

(7) 

n 

- 

1  ♦  n, 

(8) 

u 

- 

u0  +  U1 

(9) 

- 

*0  +  ♦i  • 

(10) 

Substituting  these  in  equation  (5)  we  obtain 

V2i>l  *  2n0  •  Vil; j  +  2k2n,(r)  -  0  (II) 

which  has  the  solution 

li'j(r)  - 


_ k 

2ttu 


2  fnl(r,)u0(r')  exp(ik|r  -  £'|) 

0<S>  I  Ir-r'| 


dV' 


(12) 


where  the  integration  is  carried  out  throughout  the  entire  space  occupied  by  the  medium, 
and  r'  is  a  vector  representing  a  point  in  that  space. 
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For  a  spherical  wave 


and 


(r) 


uQ(r)  -  U  exp(ikr)/r 

t2r|n,(r’)  exp[ik(r'  -  r)]  exp[ik|r  -  r*  |]  ^ 


2ir 


r'  r  -  r' 


03) 


04) 


Resolving  the  r  into  components  x  ,  along,  and  y  and  z  perpendicular  to  the  line  of 
sight,  and  neglecting  second  order  terms  in  y  and  z  ,  this  reduces  to 


4-,  (r) 


n  (r?)  exp\ik 


'(x  -  x*) 


'(x  ■  x1) 


dV1 


(15) 


Multiplying  these  equations  by  their  complex  conjugates  and  taking  the  ensemble 
average,  we  obtain  -  after  considerable  manipulation  -  the  following  relations  between  the 
two-dimensional  spectra  of  the  phase  and  amplitude  structure  functions  in  a  plane  perpen¬ 
dicular  to  the  optical  path,  at  its  far  end,  and  the  refractive  index  variations  along 
that  path. 

fa«)  -  2„k2  j  (k)2  c2(a„<°)  (!!!)  4,  .  06, 

0  v 

Fs«)  -  2„2  /  (i)2  •  <”> 

0 

) 

where  n  is  distance  measured  along  the  optical  path 
L  is  its  total  length 

is  the  value  of  <t>  for  C2  *  1  . 
n  n  n 

The  amplitude  and  phase  structure  functions  may  be  computed  from  these  expressions, 
by  means  of  equation  (4),  but  a  much  simpler  expression  is  obtained  if,  instead  of 
evaluating  these  individually ,  we  calculate  their  sum,  which  is  known  as  the  wave  struc¬ 
ture  function. 

00 

Dw(p)  «  4ir  /['  -  JoH[Fa  (K)  ♦  Fs(K)]lCdK 
0 

-  8‘2’2  /['  -  V*»]'  /<<5>  (k)2  <0>  (t)  m 

0  0 
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!0/3T(2/3)k2 
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/  (*  *■  / 


[l  -  J0(Kp)]K"^dKdn 


r(2/3)rci/6)  k2p!  I  r2 


3^?r(n/6) 


(n)dn 


2.9I4k2p^ 


(n)dn 


08) 


In  our  application,  p  is  very  small  compared  with  L  and  the  variation  in  ampli¬ 
tude  across  the  aperture  is  small  compared  with  the  variation  in  phase,  so  D  may  be 

w  j  2  ]  3 

used  in  place  of  with  little  error.  Further,  there  is  some  experimental  evidence  *  9 

that  phase  effects  saturate  less  quickly  than  those  involving  amplitude  or  intensity,  so 
the  Kolmogorov  approximation  is  more  nearly  valid. 

2.3  Modulation  transfer  function  for  long  exposures 

The  imperfections  in  the  image  produced  by  a  lens  may  be  described  in  terms  of  its 
amplitude  spread  function  which  may  be  calculated  from  the  Fourier  transform  of  the  complex 
amplitude  of  the  wave  front  at  its  exit  pupil. 

It  follows  that  the  intensity  spread  function  is  given  by  the  Fourier  transform  of 
the  autocorrelation  function  of  the  complex  amplitude  in  the  exit  pupil,  and  that  the 
autocorrelation  function  itself,  suitably  normalized,  gives  the  modulation  transfer  func¬ 
tion  of  the  lens,  see  for  example  Ref  14. 

A  parallel  result  may  be  obtained  in  the  case  of  the  imperfections  in  the  image 
formed  by  the  lens  as  a  result  of  the  random  disturbances ,  due  to  turbulence,  in  the  wave 
front  arriving  at  the  entrance  pupil  of  the  lens,  and  it  may  be  shown1 ^  that  the  modula¬ 
tion  transfer  function,  averaged  over  a  long  period  of  time,  arising  from  this  cause  is 
simply  the  spatial  coherence  function  in  the  pupil,  ie 

(M^f))  =  <E>"'  (<g(r  +  FXf)#*(r)>  .  (19) 

where  F  is  the  focal  length  of  the  lens 
1  is  the  wave  length  of  the  light 
f  is  the  spatial  frequency  in  cycles/unit  distance 
^(r)  *  expfikr  -  iwt  +  A(r,t)  +  iS(r , t)](E(r))  ^ 

hence 

<Me(£)>  «  (exp[A(r  +  FXf.t)  +  A(r)  +  iS(r  +  FXf)  -  iS(r)j) 

The  quantity  in  square  brackets  is  normally  distributed  with  mean  2<A)  ,  the  variance  of 
the  real  part  is  2a2  *  PA(F^f)}  and  the  variance  of  the  imaginary  part  is 

2 

2cjg{f  -  pg(FAf)j  where  and  are  the  amplitude  and  phase  autocorrelation  function^ 


and 


<M£(f)> 


0  . 


-  expj^2<A>  +  2o^jl  +  PA(FXf)|  -  2Cg{l  -  ps(FXf)}j2  <A>  +  2a*  = 

(M^f))  »  exp£-  2 cr*  +  2a^jl  +  pA(FAf)}  -  2cjgJl  -  p(FXf)[J  . 

By  making  use  of  the  properties  of  log-normal  random  variable  (see  eg  Ref  1 ,  sec- 


tion  6.9)  this  can  be  reduced 

to 

<  (f  )>  = 

exp 

= 

expj 

[-  jDA(FXf)  -  jDs(FXf)J 

% 

exp| 

[-  !D„(F»f>]  . 

(20) 

This  expression  represents  the  image  degradation  due  to  turbulence  when  the  exposure 
time  is  long  compared  with  period  of  the  turbulence,  and  includes  the  effects  of  both  the 
instantaneous  spread  of  the  image,  and  its  movement  in  the  focal  plane, 

(M^(f)>  may  be  expressed  in  a  dimensionless  form  as 

=  expf-  f\\ 


where  f  -  f/f„, 

t  T 

and  fT  =  0.0879(A^/  F) 

The  form  of  the  function  is  shown  in  Fig  la,  together  with  a  Gaussian  distribution 
which  has  the  same  value  at  f  *  1  for  comparison.  The  two  curves  are  superficially 
somewhat  similar  in  appearance,  but  there  is  one  important  difference:  the  second  deriva¬ 
tive  of  the  turbulence  curve  is  infinite  at  zero  frequency.  The  consequences  of  this 
will  be  mentioned  later. 

An  interesting  limiting  case  occurs  when  the  upper  limit  of  the  integral  is  extended 

2 

into  a  region  in  which  C  becomes  negligible.  Under  these  circumstances,  the  integral 
,  n  s 

becomes  proportional  to  L“3  ,  so  fT  is  proportional  to  L  for  a  given  angle  of  path, 
and  the  scale  of  the  MTF  function,  expressed  in  cycles  per  unit  length  at  the  target  is 
independent  of  L  . 

2.4  Modulation  transfer  function  for  short  exposures 

The  image  spread  averaged  over  long  periods  may  be  regarded  as  the  convolution  of 
the  instantaneous  spread  with  the  probability  distribution  of  the  position  of  its  centroid. 

Using  the  Fourier  shifting  theorem,  M  the  modulation  transfer  function  corres- 

s 

ponding  to  spread  relative  to  the  instantaneous  centroid  of  the  image,  which  is  displaced 
by  a  from  its  mean  position,  is  related  to  the  long  period  average  MTF  by 
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and  the  mean  values  by 


M  (f)  =  M  (f)  exp[-  2irict  .  f] 


(M  (f))  =  (M  (f))(exp(-  2iria  .  f)) 


=  (Mg(f))  exp[2ir(a  .  f>3 
=  (Ms(£))  exp[-  ir2f2(a  .  a)] 

cx  is  given  by  the  mean  rate  of  change  of  phase  over  the  pupil,  multiplied  by  F/k  »  so  in 
the  one  dimensional  case 

<£  «  £i>  =  (F/pk)2Dg(p) 

where  p  is  the  width  of  the  pupil. 

In  two  dimensions,  we  would  expect  the  value  to  be  approximately  twice  this,  and  it  can  be 
shown*  that  for  a  clear  circular  aperture  of  diameter  d 


4(x)f  W  s 

*  °-975  (0Dsw 


<  a  .  a)  =  2  * 


The  correcting  factor  is  normally  taken  to  be  unity  so 


\  ( f  ))  =  (M^Cf))  exp£~  ^2(F/dk)  ^Dg(d)^ 


re-writing  in  terms  of  the  cut  off  frequency  of  the  lens 


kd.^irF 


and  since 


<iyn)  =  <Mg(n)  exp£-  Hf/f0)2Dg(d)] 

-  <Ms(f)>  exp i(f/f0)^Ds(FXf)] 


°S  *  Dw 


(Ms(f)>  *  exp£-  i|l  -  (f/f0)^Dw(FXf)J  . 


(f )  ,  unlike  Mf(f)  depends  on  the  aperture  of  the  lens  being  used  to  collect  the  radia¬ 
tion,  so  it  cannot  strictly  ve  regarded  is  a  transfer  function  for  the  optical  path,  but 
only  as  a  multiplying  factor  ^e  p  lied  to  the  MTF  of  the  lens.  It  has  no  significance 
for  values  of  f  >  where  the  MTF  of  the  lens  is  identically  zero. 
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Equation  (21)  may  be  expressed  in  the  form 


(Mg(f)>  =  exp{-  57 .52t(i  - 


where  f  =  t/t 
r  0 


and  T  is  a  dimensionless  parameter 


T  =  (d*/X2)  J  (jij  C2(n)dn 


all  the  lengths  being  expressed  in  the  same  units. 

The  minimum  value  of  <M  (f)  >  occurs  at 

s 


«.  ■  (!) 


*  0.5787 


and  has  the  value 


M  .  =  exp(-  3.853T) 

min 


This  is  probably  the  most  convenient  single  parameter  for  expressing  the  effect  of  turbu¬ 
lence  in  short  exposures. 

The  form  of  <Ms(f)>  for  a  range  of  values  of  T  is  shown  in  Fig  2  together  with 
the  corresponding  values  of  (M  (£))  for  comparison. 

2 . 5  Intermediate  exposures 

The  way  in  which  the  form  of  the  turbulent  "modulation  transfer  function"  varies 
l  exposure  time  depends  on  the  specific  conditions.  Hufnagel *  states  that  for  photo¬ 
graphic  systems  with  exposure  times  longer  than  several  milliseconds  image  motion  contri¬ 
butes  significantly  to  the  overall  blur.  On  the  other  hand,  Roddier  and  Roddier^,  who 
were  concerned  primarily  with  astronomical  images,  found  that  the  difference  between  the 
image  power  spectra  at  zero  and  0.02  second  exposure  time  is  negligible,  and  even  at 
0.1  second  exposure  the  difference  is  small  except  close  to  the  cut  off  frequency. 

Further  work  is  required  on  this  problem,  but  as  the  exposure  times  used  in  air  to  ground 
photography  are  of  the  order  of  a  very  few  milliseconds,  the  short  exposure  form,  equation 
(21)  may  usually  be  employed. 

2 . 6  Line  spread  function 

A  line  spread  function  representing  the  effect  of  turbulence  in  long  exposures  may 
be  obtained  by  Fourier  transforming  equation  (20),  and  is  shown  in  Fig  lb  as  a  function  of 


*  f  x 
T 


It  falls  off  very  slowly  at  large  distances,  and  because  the  function  transformed  has  an 
infinite  second  derivative  at  zero  frequency,  its  standard  deviation  is  infinite.  If  a 
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characteristic  dimension  is  required,  the  quartile 

-  0. 1 50/ f 

q  T 

or  the  standard  deviation  of  the  equivalent  Guassian  distribution 

0G  =  l/(/2irfT)  ^  0.225/fT 

may  be  used.  Other  linear  parameters  may  be  obtained  from  Fig  lb.  The  effects  of  turbu¬ 
lence  in  short  exposures  cannot  be  expressed  directly  in  terms  of  a  line  spread  function 
because  the  Fourier  transform  of  equation  (21)  does  not  converge,  but  the  combined  line 
spread  function  of  a  turbulent  path,  and  the  receiving  lens,  may  be  computed  by  Fourier 
transforming  the  product  of  their  modulation  transfer  functions. 

The  combined  MTFs  of  a  diffraction  limited  lens,  and  a  number  of  turbulent  optical 
paths  are  shown  in  Fig  3,  and  the  corresponding  line  spread  functions  are  plotted  against 

x  *  f~x 
r  0 

in  Fig  4.  The  corresponding  long  exposure  cases  are  also  shown  for  comparison. 

3  ESTIMATING  THE  STRUCTURE  CONSTANT  OF  REFRACTIVE  INDEX 

3 . 1  General  approach 

2 

The  way  in  which  C~  varies  with  place  and  time  has  been  measured  by  a  number  of 

1-9  .  n  , 

workers  .  As  might  perhaps  be  expected,  the  variations  are  large  and  erratic,  and 

there  would  seem  no  prospect  of  predicting  them  in  detail  for  some  point  in  space  on  a 

given  occasion. 

For  the  purposes  of  estimating  the  performance  of  reconnaissance  systems,  and  of 
optimizing  this,  however,  all  we  require  are  weighted  average  values  along  an  optical  path 
for  a  given  time  of  day,  season,  altitude  and  broadly  specified  geographical  area  and 

weather  situation.  It  may,  in  time,  be  possible  to  accumulate  sufficient  statistical  data 

2  .  .... 
on  to  estimate  these  directly,  but  in  view  of  the  difficulty  of  the  measurements  this 

seems  a  rather  remote  possibility,  and  for  the  present  we  will  have  to  estimate  the  mean 

2 

values  of  indirectly  from  meteorological  parameters  for  which  detailed  statistics 

are  generally  available. 

Such  estimates  are  unlikely  to  be  very  accurate,  but  the  large  random  errors  which 
will  occur  will  probably  swamp  any  bias  introduced,  and  they  should  be  adequate  for  broad 
system  design  work. 

The  dominant  meteorological  factors  involved  differ  at  high  and  low  altitudes  and 
the  two  cases  will  be  discussed  separately. 

3. 2  Lower  atmosphere 

The  local  variations  in  the  refractive  index  of  the  atmosphere  in  turbulence  arise 
primarily  from  the  variations  in  temperature,  and  are  related  to  them  by  the  equation 
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4  =  78  x  io“6  4  >  {22) 

dT  2 

where  P  is  the  local  pressure  in  mb 
./■  is  the  absolute  temperature. 

For  rough  statistical  work,  values  of  P  and  •'f  based  on  the  ICAO  standard 
atmosphere^  may  be  used. 

In  the  lowest  layers  of  the  atmosphere  turbulence  is  primarily  due  to  insolation. 
It  is  therefore  convenient,  first  to  estimate  the  structure  constant  of  temperature,  and 
then  to  derive  the  structure  constant  of  refractive  index  from  it  using  equation  (22) 
and  the  relation 


c2  =  (dnf  C2 
n  \df)  .i 


A  number  of  formulae  have  been  given  for  the  variation  of  Cj  with  altitude  in  the 
boundary  layer.  The  simplest  of  these  is  a  negative  four  thirds  power  law^*  ^ 

C  r  «  h  3  .  (24) 

However,  this  must  break  down  at  very  low  altitudes  as  not  only  does  it  give  infinite 

2  2  2 

values  of  and  at  zero  height,  but  the  weighted  integral  of  from  zero  to  any 

finite  height  is  also  infinite. 

|9  .  .  .  .  A 

Wyngaard  ct  ai  have  given  a  modified  formula  which  varies  as  h  at  very  low 

height,  but  reverts  to  h  3  further  up,  the  transition  taking  place  at  a  height  of  a  few 

metres.  This  expression  is  also  infinite  at  zero  height,  but  the  integral  remains  finite, 

and  tends  to  zero  as  the  path  of  integration  tends  to  zero. 


A  convenient  form  of  the  formula  has  been  given  by  Hufnagel  . 


r(h)  - 


2  x  10  Q’h  3 

3*|$ 

1 400u^ 

'  +  hQ 


where  Q  is  the  upward  convective  heat  flux  W  m 
h  is  the  height  in  metres 

u^  is  a  characteristic  frictional  velocity  ms1 

2 

If  Q  and  u  can  be  estimated,  C  can  be  evaluated  for  each  point  on  the 
*  n 

optical  path  using  equations  (23)  and  (25),  and  1)  ,  the  wave  structure  function  at  the 

w 

end  of  the  path  computed  by  the  integrating  of  equation  (18).  This  integration  must  be 
performed  for  each  Q  ,  u*  and  path  configuration. 

Q  is  dependent  on  such  factors  as  solar  elevation,  cloud  cover  and  the  type  of 
terrain.  As  a  rough  guide1  we  may  take 

Q  *  sin  -  50  (26) 
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where  £  is  the  solar  elevation  and  typical  values  of  are  as  shown  in  Table  1 . 

When  the  right  hand  side  of  equation  (26)  is  negative,  Q  may  be  taken  as  zero. 


Table  1 

2 

Qq  in  watts  per  metre 


^  '~--^Sky  condition 

Type  of  terrain^~^^__ 

Continuously 

clear 

Overcast 

Dry  sand  or  lava 

500 

200 

Dry  fields  of  brush 

400 

150 

Wet  fields 

200 

70 

For  intermediate  amounts  of  cloud  we  may  estimate  by  linear  interpolation. 

The  use  of  the  words  ’continuously  clear*  in  Table  1  should  be  noted,  as  changes  in  tur¬ 
bulence  always  lag  somewhat  behind  the  meteorological  factors  producing  them.  A  similar 
slight  asymmetry  will  occur  between  morning  and  evening  conditions. 

u*  is  typically  an  order  of  magnitude  smaller  than  the  local  wind  speed.  Its 
2  .  .  . 

effect  on  is  only  significant  at  low  heights,  and  its  contribution  to  the  weighted 

integral  even  smaller,  if  the  detector  is  well  above  the  transition  point  as  will 

normally  be  the  case  in  air  to  ground  reconnaissance.  For  most  of  our  applications  we 

3 

have  therefore  used  a  mean  value  of  U*/Q  which  gives  a  break  point  between  the  two  forms 
at  a  height  of  a  few  metres,  say  three,  and  written 


C^h)  = 


2  x  io 


;  *  If 


This  has  the  great  advantage  that,  the  integral  in  equation  (18)  can  be  expressed  as  a 
£ 

product  of  Q1  and  a  function  of  height  which  may  be  pre-computed  and  fitted  by  an 
empirical  function  for  routine  use  as  will  be  described  in  section  4.2. 

3. 3  Upper  atmosphere 

2 

Above  the  boundary  layer,  the  mean  value  of  falls  off  rather  more  slowly  than 

at  low  levels,  though  a  sharp  local  maximum  occurs  in  the  neighbourhood  of  the  tropopause, 
which  is  often  a  region  of  high  wind  shear. 

I  20 

Hufnagel  *  has  given  an  empirical  formula  expressing  the  structure  constant  of 
refractive  index  in  the  region  above  the  first  strong  temperature  inversion  as  a  function 
of  height  and  of  the  root  mean  square  wind  speed  averaged  over  the  height  band  5-20  km, 

C*(h)  =  8.2  x  ]0“56U2h'°exp(-h/1000)  +  2.7  x  10~16  exp(- h/l500)m”l  (28) 


where  h  is  in  metres 

U  is  the  average  rms  wind  speed  in  the  height  band  5-20  km  in  metres  per  second. 
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Values  of  U  for  any  given  area  and  season  may  be  deduced  from  the  maps  in  Upper 
Winds  over  the  World"  .  Hufnagel  suggests  27  m  s  as  typical,  but  this  appears  to  have 
been  based  on  experiments  in  an  area  of  higher  than  average  upper  winds. 

"“17  — ^ 

For  U  =  27  ,  equation  (28)  gives  a  peak  value  at  h  *  1 1  km  of  about  10  m  ^  . 
This  is  numerically  quite  close  to  the  peak  value  shown  in  Fig  5  of  Hufnagel* s  earlier 
paper  and  Fig  6  of  Hufnagel  and  Stanley's  paper  ,  but  as  these  are  expressed  in  units 

of  cm  they  represent  a  much  larger  turbulence.  However,  as  far  as  can  be  deduced  from 
the  small  scale  graphs,  the  peak  is  also  very  much  narrower  than  that  given  by  equation 
(28),  and  it  may  be  that  the  explanation  of  this  discrepancy  is  that  the  earlier  graph 
represents  a  typical  profile,  while  equation  (28)  represents  an  ensemble  average  of  all 
occasions,  and  the  peak  is  therefore  flattened  by  the  variations  in  the  height  of  the 
tropopause. 

3.4  Intermediate  heights 

The  form  of  the  transition  between  the  two  types  of  turbulence  structure  is  highly 
dependent  on  the  specific  weather  situation,  but  unless  it  is  possible  to  forecast  the 
thickness  of  the  boundary  layer,  the  best  that  can  be  done  is  to  flare  the  curves  into 
each  other  in  the  transition  region,  so  as  to  approximate  to  which  ever  is  the  greater. 

The  simple  form 

CN  =  (Cn)  CXp(~  °-00001h2>  +  (Cn)  (29) 

'  '  'low  \  /high 

fits  the  data  with  a  peak  error  of  less  than  5%  for  a  Q  of  100.  There  is  no  need  to 
attenuate  the  |cfT  j  term  at  low  altitudes  as  it  becomes  small  compared  with  |C^  ) 

\  N/high  V  N/lov 

unless  both  are  negligible. 

2 

Cn(H)  may  be  expressed  as 

C^(H)  =  Cj (H)  +  C2(H)U2  +  C3(H)Q3 

and  the  forms  of  the  coefficients  C|,  and  are  shown  in  Fig  5. 

4  COMPUTATIONAL  FORMULAE 
4 . 1  Modulation  transfer  function 

The  expressions  for  long  and  short  exposure  modulation  transfer  function  can  be 
combined  in  the  form 

{^f\ 

L 

where  IL  -  J  (tfc  2(n)dn 
0 


a  i 
X  3 (Ff )3 I 


s  is  0  for  slow  shutter  speeds 
]  for  fast  shutter  speeds. 


065 


15 


Rewritten  in  terms  of  more  convenient  units  we  have 


[-  57528{<  "  0 • 0 1 S Kf«)1  U(Fmmf,nm)5 


M(f  )  *  exp 

mm 


where  N  is  the  f /number  of  the  lens 

F  is  the  focal  length  of  the  lens  in  millimetres 

mm 

X  is  the  wave  length  of  the  light  in  nanometres 
nro 

f  is  spatial  frequency  in  cycles  per  millimetre 
mm 


or  in  angular  terms 


M(w  )  =  exp|-  0.5753  x  10 

mr  | 


v  1 

„  /x  a)  v  M 

10  i-o-i.lW  to 

\  mm  / 


where  oj  is  the  spatial  frequency  in  cycles  per  milliradian 
mr 

d  is  the  diameter  of  the  lens  in  millimetres, 
mm 

The  latter  form  is  used  in  the  subroutine  to  be  described,  as  it  involves  one  fewer 


variable . 

4,2  Evaluation  of  the  integral 

Consider  an  optical  path  from  a  target  at  height  to  a  camera  at  height  Hj  at 

a  plan  range  of  Rp  (see  Fig  6) .  Then  if  we  assume  that  over  this  optical  path  Cn 
may  be  treated  as  a  function  of  height  only, 


(if  <  (:>d" 


_ _ 1 

•,  -  v2  *  rp  *  t(»,  -  vW 


where  £  =  n/L  . 

For  paths  with  a  small  height  range  the  integral  may  be  approximated  by  one  of  the  follow- 
ing  generalized  Gauss  quadrature  formulae 


[0. 1 I95C^(0.574Hq  +  0.4261^) 
0. 2555C2 (0. 1 32Hq  +  0.868Hj) 


which  are  exact  if  C2  varies  linearly,  or  as  a  cubic  in  H  respectively,  over  the 
n 

required  interval. 
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For  wider  ranges  of  H  the  weighted  integral  becomes  complex  in  form  and  is  best 
evaluated  piecewise  using  variable  step  length. 

For  paths  starting  at  the  ground  we  may  write 


1  + 


(Rp/H,r 


where  I 


H 


/ 

0 


(h)dh 


This  special  case  is  of  very  frequent  occurrence,  and  to  save  time  in  routine  analy¬ 
sis  the  integral  I  has  been  pre-coxnputed  for  1  <  H  <  10^  metres,  and  fitted  by  an 
empirical  function. 

d 

IH  =  I,  (Hj)  +  iti2(h,)  +  q3i3(h,)  . 

Simple  polynomial  fits  to  1^,  1^  and  1^  were  not  successful,  but  it  was  found 

that  their  logarithms  could  be  expressed  as  the  difference  of  a  polynomial,  and  a  root  of 

a  polynomial  in  log  Hj  .  The  coefficients  of  the  two  polynomials  for  each  function  were 

determined  by  the  iterative  use  of  the  ICL  curve  fitting  sub-routine  F4CF0RPL  and  are 

shown  in  Table  2,  together  with  their  mean  square  error  weighted  to  take  account  of  the 

relative  importance  of  the  three  terms  at  any  height.  Two  versions  of  I0  are  given. 

2  .  1 

In  view  of  the  probable  error  in  estimating  U  ,  the  simpler  form  will  often  be  adequate. 


Table  2 


Empirical  expressions  for  1^,  and 
10  +/-  9.293  -  0. 3977y  -  0.0l883y2  -  Vo.8 1 51  -  I.805y  +  I.802y2  )  (2. 

3  -  J60.35  -  108. Oy  +  59.30y2) 


.5  Z) 


10  +  ^-  13. 


04  +  2.886y  +  0.2075y  +  0.l520yJ 


I2  -  ioH 


I-  13.53  +  3.453y  +  0.2414y2  +  0.1045y3  +  0.007853y 


*Jrj\2  -  10009y  +  1 5544y2  -  1 1 629y3  +  3527y4 


35  -  108. Oy  +  59. 30y  ) (8Z) 

A\ 

(2%) 


0+^-  14. 


20  -  0.3068Y 


-  0.1317Y2  -  J 1 .445  -  I.177Y  +  0.2405Y2) 


(3%) 


where  y 
Y 


log  (Hj / 1000) 
log  (Hj/10) 


H 


height  of  sensor  in  metres 
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4.3  FORTRAN  subroutines 

The  equations  described  in  the  previous  section  have  been  implemented  in  four  sub¬ 
routines  and  an  auxiliary  function  segment  written  in  extended  XCL  1900  FORTRAN,  which  are 

24 

included  in  the  * PHOTRAN*  subroutine  library 

They  are  as  follows: 

SUBROUTINE  TIH(H1 ,Q,U,TI) 
which  evaluates  I 

H 

SUBROUTINE  TIL(H0,H1 ,RP,Q,U,TI) 
which  evaluates  1^ 

SUBROUTINE  CSN(H,Q,U,C) 
which  evaluates 

n 

SUBROUTINE  MTFT (NEW,FA, S,D,WL,H0,H1 ,RP,Q,U,T) 
which  computes  the  modulation  transfer  function. 

FUNCTION  CSNFUN(RSF) 

5  2 

which  provides  values  of  n5  (n)  for  evaluating  IT  for  paths  which  do  not  start  at 

n  Li 

ground  level. 

The  variables  used  in  the  CALL  statements  are  as  follows 

HO  target  height  in  metres 

HI  camera  height  in  metres 

RP  plan  range  in  metres 

WL  wave  length  in  nanometres 
D  lens  diameter  in  millimetres 

S  shutter  speed  parameter,  0  for  long  exposure 

1  for  short  exposures. 

The  subroutine  will  accept  intermediate  values. 

NEW  is  set  to  1  when  a  new  optical  path  is  required,  and  to  0  when  only  the  fre¬ 
quency  is  changed, 

-2 

Q  upward  convective  heat  flux  in  watts  metre 

U  the  average  rms  wind  speed  in  the  height  band  5-20  km,  in  metres  per  second 
FA  spatial  frequency  in  cycles  per  milliradian. 

If  FA  is  set  to  -I  ,  the  frequency  at  which  the  short  exposure  'MTF*  is  least 
will  be  substituted. 

T  Modulation  transfer  function  at  a  frequency  FA 
RSF  -  n/L 

HO,  Hi,  Q  and  U  are  passed  from  TIL  to  CSNFUN  via  a  common  block 
C0MM0N/TURBDATA/H0,H1,Q,U  . 

Detailed  listings  are  given  in  the  Appendix. 
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One  point  that  should  be  noted  is  that  if  subroutine  TIL  is  called  with  HI  <  1 
metre  it  will  be  set  to  1  metre  and  this  value  will  be  used  in  the  computation,  and 
returned  to  the  calling  segment  at  exit* 

5  DISCUSSION 

The  theoretical  difficulties  and  the  scarcity  of  experimental  data  have  been 
stressed  throughout  this  Report,  and  there  seems  no  immediate  prospect  of  predicting  the 
effect  of  turbulence  on  the  image  quality  of  an  individual  exposure. 

The  methods  presented  here,  however,  should  provide  a  convenient,  and  adequately 
accurate  way  of  introducing  the  average  effect  of  turbulence,  in  some  area  and  broadly 
defines  meteorological  situation,  into  mathematical  models  for  use  in  the  analysis  and 
optimization  of  reconnaissance  systems,  both  photographic  and  electro-optical.  Some 
detailed  examples  of  their  application  will  be  given  elsewhere. 

Perhaps  the  most  important  result  to  emerge  is  that  there  is  a  large  difference  in 
the  performance  obtainable  in  short  and  long  exposures.  High  aperture  optical  systems  are 
to  be  preferred,  not  only  because  they  permit  shorter  exposures,  and  hence  reduce  or 
eliminate  the  turbulent  image  movement,  but  also  because  the  residual  static  blurr  is 
aperture  dependent. 

The  exact  form  of  transition  between  the  long  and  short  exposure  cases  is  one 
field  which  would  repay  further  investigation. 
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Appendix 

LISTINGS  OF  SUBROUTINES 

SURROUT  IMF  C5N<H,Q,U,C) 

T*?1A- AS 

P=?2A. ?4*FXP<  0001  5765*  ni000-H)1 

!f<h.gt.  hoooigoto  i 

T*288.18-  0065*H 

P*1 01 3. 2*PXP(5. 25766*4100(1-. 0000225 5 5 35 *H) ) 

1  G*<  .  f*00O79*P/ ^t*t)  )  **2 
C  =  2.7f-16*EXP(-H/1500) 

1  +U*U*R. ?F-54*M**1O*FXO(-.001*M) 

1  +G*EXP<-.00001*H*h)*.00?*<0/H>**(4./3.1/(1.*3./H)**<2./3.) 

RETURN 
FND 


SUROOUMNF  t  ih  (HI  .0  »U  ,TJ  ) 

DAT*  ag . »1 , A? /-« . 20 3. - .39  77, - . 01 883  / 

DATA  80.8* .8?/.  8151  ,-1.805*1 .80?/ 

DATA  CO.Ci ,C?.r3,C4/-1 3.53, 3. 453.  2414, . 1 045 , . 007853 / 
DATA  DO."- ,P?,D3,D4/ ?71 2.-1 0009, 15544.-1 1029,3577/ 
DATA  FQ.F1 , E?/-1 4 . 2 3008#- . 1 31 7/ 

DATA  FO. FI , F2/1 . 445.-1 .177, . 2405/ 

2*41001 o (Hi ) -1 . 

Vs 7-2. 

Tt*10**(«n+Y*<Al*V*A2)-5ORT(P0+V*(A1+Y*82))) 

1+U*U*1f'**<CO-*’V*(C1+Y*(C?*Y*(C3*V*e4))> 

1-  (00*  V*  1 0l  +  V*(r>2*V*  <p3*Y*D4)  )))**(.  25)  > 

1*0**(4./3.)*10**<EO'»2*<FlT7*F?)-FOPT(FO  +  Z*(Fl*7*F2))) 
RETURN 
END 


SURROUTINF  MTFT (NEW  ,FA  ,5  ,D  ,WL  ,H0  ,H1  .PP.O.IJ.T) 
IF  (NFW.  FQ.  1  )CA(Ll  T t  L ( HO . Hi , RP . 0 , U , T t ) 

TF(FA . IE. -0.000001 >FA*578.7*D/WI 

TsFXO(-.5753F10*(1-.1*S*(WL*FA/D)**<1./3)1 
1  *Wl**'-1./3.>*FA**<5./3.)*Tl) 

RETURN 

END 
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SUBROUTINE  TlltHO.HI ,RP.O,U,Tl) 

COMMON  /  TURBO  AT  A /HHO  *HHl  ,OQ,UU 

external  crnfijn 
IF (HI . IT.l .  )M1*1 . 

HHO>HO 
HH1 sHI 
0Q=0 
IJU»U 

PS  =  SOOTf(n1-HO)**2'»RP*PP) 

I F (HO . LT . 0 . 001 ) GOTO  300 

CALI.  CSN(.574*hO+.A?6*h1*Q»U,CO) 

CALI  CRN ( . 1 32*H0+. A6B*H1 ,Q,U.C1 ) 

CBR,CB=.11O5*C0+.2S5S*C1 

IF(ABS(H1/HO-i.).6t.0.3)CAll  F4 1  NtSMP ( 0 . ,1 . , CSNCUN. . 001 *CB » 3 , CBS > 

TI=CRB*RS 

RETURN 

3()0  CALL  t  |  H  f  N1  ,Q,U»TH1 
T I sTH*RS  /Hi 
RETURN 
END 

FUNCTION  CSNFIINLRSF) 

common/tmrboata/ho.hi ,0,U 

HaHO*(H1-HO)*PRF 
CALL  CSN(H, 0,11,0 

CSNFUN*PRF**(S. /3. W  \ 

return 

END 
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